Geospatial Analysis
Davit Kajaia
1 Part 1
1.1 Libraries
library(tidyverse)
library(ggspatial)
library(maptools)library(dplyr)
library(ggthemes)
library(spatstat)
library(classInt)
library(readxl)
library(knitr)
library(sf)
1.2 Loading and tidying data
setwd("C:/Users/User/Documents/OneDrive/Desktop/R/spatial")ge_admin_url <- "https://data.humdata.org/dataset/3ee95199-2dfe-40fc-b9bd-b44cc8c91024/resource/ac2b2747-18a2-43c9-abbb-4ab2c5c539a3/download/geo_adm_geostat_20191018_shp.zip"
download.file(ge_admin_url, destfile="ge_admin.zip")
unzip("ge_admin.zip")
file.remove("ge_admin.zip")ge_adm_0 <- st_read("geo_admbnda_adm0_geostat_20191018.shp")## Reading layer `geo_admbnda_adm0_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm0_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 13 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 40.00834 ymin: 41.05455 xmax: 46.72672 ymax: 43.5865
## geographic CRS: WGS 84ge_adm_1 <- st_read("geo_admbnda_adm1_geostat_20191018.shp")## Reading layer `geo_admbnda_adm1_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm1_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 13 features and 16 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 40.00834 ymin: 41.05455 xmax: 46.72672 ymax: 43.5865
## geographic CRS: WGS 84ge_adm_2 <- st_read("geo_admbnda_adm2_geostat_20191018.shp")## Reading layer `geo_admbnda_adm2_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm2_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 70 features and 19 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 41.54717 ymin: 41.05455 xmax: 46.72672 ymax: 43.25688
## geographic CRS: WGS 84ge_pop<-read_xlsx("geo_admpop.xlsx",sheet=3)## New names:
## * `` -> ...8ge_pop_2 <- ge_pop%>%
select(ADM2_PCODE, POP_Total,F_Total,M_Total)
ge_adm_2_ii <- ge_adm_2 %>%
select(ADM2_PCODE, ADM2_EN)
ge_adm_2_pop <- left_join(ge_adm_2_ii, ge_pop_2, by="ADM2_PCODE")1.3 Adding data for Tbilisi
ge_adm_1_TBI <- ge_adm_1 %>%
select(ADM1_PCODE, ADM1_EN) %>%
rename(ADM_PCODE = ADM1_PCODE,
ADM_EN = ADM1_EN) %>%
filter(ADM_EN == "Tbilisi")
ge_adm_1_TBI <- ge_adm_1_TBI %>%
mutate(F_Total=634743,M_Total=536357,POP_Total=sum(634743,536357))
ge_adm_2_pop <- ge_adm_2_pop %>%
rename(ADM_PCODE = ADM2_PCODE,
ADM_EN = ADM2_EN)
ge_adm_2_pop <- rbind(ge_adm_2_pop, ge_adm_1_TBI)1.5 Plotting
ggplot()+
theme_minimal(base_size=8.5)+
theme(axis.title = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size = 13,face="bold",hjust=0.5))+
geom_sf(data = ge_adm_0, fill = "grey", colour = "grey60", size = 0.1)+
geom_sf(data = ge_adm_2_pop, aes(fill = Share*100),colour = "grey60", size= 0.1)+
scale_fill_gradientn(colors=c("#004648","#E6E2BC","red")) +
labs(title = "Gender Balance in Georgia",
fill = "Share of female \npopulation,%",
caption = "data: The Humanitarian Data Exchange \nvisual: D.Kajaia")+
annotation_scale(location = "bl", style = "ticks") +
annotation_north_arrow(location = "tr", width = unit(0.5, "cm"))
ggsave("Gender_balance_in_Georgia.png", dpi = 600, height = 4, width = 7, units = "in")1.6 Conclusion
ge_adm_2_pop %>% summarize(sd=sd(Share),mean=mean(Share))## Simple feature collection with 1 feature and 2 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 41.54717 ymin: 41.05455 xmax: 46.72672 ymax: 43.25688
## geographic CRS: WGS 84
## sd mean geometry
## 1 0.01318941 0.5123185 POLYGON ((46.39876 41.47027...From the map we can see that the sex ration is not equal - it is more female-biased. Main reason for that can be that females live longer than males.
As we can see there is some difference in gender share between districts, but Standard deviation of gender share is very small - only 1 %, therefore mostly these differences are uxplainable. It seems that female share is higher in some cities with higher population (for example in Tbilisi). This may be the case due to selective abortion. Selective abortion should be lower in big cities because people in big cities should be more educated.
2 Part 2
2.1 Libraries
library(tidyverse)
library(ggspatial)
library(readxl)
library(knitr)
library(curl)
library(sf)
2.2 Importing data
setwd("C:/Users/User/Documents/OneDrive/Desktop/R/spatial")
url <- "http://aasa.ut.ee/Rspatial/data/FarmedAnimalsByLocation_31102018.xlsx"
destfile <- "FarmedAnimalsByLocation_31102018.xlsx"
curl_download(url, destfile)
agrAnimal <- read_excel(destfile)download.file("https://geoportaal.maaamet.ee/docs/haldus_asustus/omavalitsus_shp.zip", destfile="omavalitsus_shp.zip")
unzip("omavalitsus_shp.zip")
file.remove("omavalitsus_shp.zip")list.files(pattern = ".shp")## [1] "geo_admbnda_adm0_geostat_20191018.shp"
## [2] "geo_admbnda_adm0_geostat_20191018.shp.xml"
## [3] "geo_admbnda_adm1_geostat_20191018.shp"
## [4] "geo_admbnda_adm1_geostat_20191018.shp.xml"
## [5] "geo_admbnda_adm2_geostat_20191018.shp"
## [6] "geo_admbnda_adm2_geostat_20191018.shp.xml"
## [7] "geo_admbndl_admALL_geostat_itos_20191018.shp"
## [8] "geo_admbndl_admALL_geostat_itos_20191018.shp.xml"
## [9] "geo_admbndp_admALL_geostat_itos_20191018.shp"
## [10] "geo_admbndp_admALL_geostat_itos_20191018.shp.xml"
## [11] "gis_osm_buildings_a_free_1.shp"
## [12] "gis_osm_places_free_1.shp"
## [13] "gis_osm_roads_free_1.shp"
## [14] "gis_osm_waterways_free_1.shp"
## [15] "gps_us.shp"
## [16] "omavalitsus_20200504.shp"municip <- st_read("omavalitsus_20200504.shp") ## Reading layer `omavalitsus_20200504' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\omavalitsus_20200504.shp' using driver `ESRI Shapefile'
## Simple feature collection with 79 features and 5 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 369032.1 ymin: 6377141 xmax: 739152.8 ymax: 6634019
## projected CRS: Estonian Coordinate System of 19972.3 Cleaning data
agrAnimal$county <- NULL
agrAnimal$`admin unit` <- NULL
agrAnimal <- agrAnimal %>% rename(action_place = `action place`,
estonian_holstein_cattle = `estonian holstein cattle`,
estonian_red_cattle = `estonian red cattle`,
estonian_native_cattle = `estonian native cattle`,
beef_cattle = `beef cattle`,
y = `X koordinaat`,
x = `Y koordinaat`)
agrAnimal <- agrAnimal %>% mutate(x = as.numeric(x), y = as.numeric(y))agrAnimal_2 <- gather(agrAnimal, "key", "value", 2:8)
agrAnimal_2 <- agrAnimal_2 %>% filter(value > 0)
agrAnimal_2 <- agrAnimal_2 %>%
mutate(municip_key = str_to_lower(municipality))
municip <- municip %>%
mutate(municip_key = str_to_lower(ONIMI))2.4 Spatial join
agrAnimal_2_sf <- st_as_sf(agrAnimal_2, coords = c("x", "y"), crs = 3301)
agrAnimal_2_sf_municip <- st_join(st_transform(agrAnimal_2_sf, 3301), st_transform(municip, 3301), join = st_intersects)
st_geometry(agrAnimal_2_sf_municip) <- NULL
agrAnimal_2_sf_municip_aggr <- agrAnimal_2_sf_municip %>%
group_by(OKOOD, key) %>%
summarise(sum = sum(value)) %>%
ungroup()
agrAnimal_2_sf_municip_aggr <- agrAnimal_2_sf_municip_aggr %>%
spread(key, sum)
agrAnimal_2_sf_municip_aggr <- left_join(municip, agrAnimal_2_sf_municip_aggr, by="OKOOD")2.5 Pigs distribution in municiaplaities
ggplot()+
theme_grey(base_size=8.5)+
theme(plot.title = element_text(size = 13,face="bold"))+
geom_sf(data = agrAnimal_2_sf_municip_aggr, aes(fill=pigs), size=0.25, colour = "grey70")+
scale_fill_gradientn(colours = c("#B0B579","#5E623F"), na.value = "#B6D3C1")+
labs(fill = "N",
title = "Pigs in Estonian Municipalities",
subtitle = "Agricultural Registers and Information Board",
caption = "Author: D.kajaia")
2.6 Density
agrAnimal_2_sf_municip_aggr$area <- st_area(agrAnimal_2_sf_municip_aggr)
agrAnimal_2_sf_municip_aggr <- agrAnimal_2_sf_municip_aggr %>%
mutate(area = as.numeric(area) / 1000000) ggplot()+
theme_grey(base_size=8.5)+
theme(plot.title = element_text(size = 13,face="bold"))+
geom_sf(data = agrAnimal_2_sf_municip_aggr, aes(fill = pigs/ area), size=0.25, colour = "grey70")+
scale_fill_gradientn(colours = c("#B0B579","#5E623F"), na.value = "#B6D3C1")+
labs(fill = "N per km2",
title = "Pigs Density in Estonian Municipalities",
subtitle = "Agricultural Registers and Information Board",
caption = "Author: D. Kajaia")+
annotation_scale(location = "bl", style = "ticks") +
annotation_north_arrow(location = "tr", width = unit(0.5, "cm"))
2.7 Finding municipality
agrAnimal_2_sf_municip_aggr%>%
filter(pigs==max(pigs,na.rm=TRUE))%>%
select(ONIMI,OKOOD,pigs)## Simple feature collection with 1 feature and 3 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 576487.5 ymin: 6437907 xmax: 628026 ymax: 6497708
## projected CRS: Estonian Coordinate System of 1997
## ONIMI OKOOD pigs geometry
## 1 Viljandi vald 0899 78237 MULTIPOLYGON (((621049.2 64...agrAnimal_2_sf_municip_aggr%>%
mutate(density=pigs/area)%>%
filter(density==max(pigs/area,na.rm=TRUE))%>%
select(ONIMI,OKOOD,density)## Simple feature collection with 1 feature and 3 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 576487.5 ymin: 6437907 xmax: 628026 ymax: 6497708
## projected CRS: Estonian Coordinate System of 1997
## ONIMI OKOOD density geometry
## 1 Viljandi vald 0899 57.03777 MULTIPOLYGON (((621049.2 64...Both – the number and the density of pigs are highest in Viljandi.
3 Part 3
3.1 Libraries
library(rnaturalearth)
library(lubridate)
library(ggspatial)
library(tidyverse)
library(classInt)
library(stringr)
library(readxl)
library(knitr)
library(units)
library(sf)
3.2 Importing data
setwd("C:/Users/User/Documents/OneDrive/Desktop/R/spatial")dt<-read.csv("https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv")
countries50 <- ne_download(scale = 50, type = 'countries', category = 'cultural', returnclass = "sf")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\User\AppData\Local\Temp\Rtmp0YCM6x", layer: "ne_50m_admin_0_countries"
## with 241 features
## It has 94 fields
## Integer64 fields read as strings: POP_EST NE_ID3.3 Cleaning data
recovered <- gather(dt, "key", "value", 5:ncol(dt))
recovered$key<-str_replace_all(recovered$key,'X',"")
recovered<-recovered %>% mutate(place=paste0(Country.Region, ", ", Province.State),date=mdy(key))
recovered$place<-gsub(", $", "",recovered$place)
recovered_sf <- st_as_sf(recovered, coords = c("Long", "Lat"), crs = 4326)recovered_sf_latest <- recovered_sf %>%
filter(date == max(date)-days(1) | date == max(date)- days(8))
recovered_crd_latest <- st_coordinates(recovered_sf_latest) %>%
as_tibble()
recovered_sf_latest_2 <- bind_cols(recovered_sf_latest, recovered_crd_latest)
recovered_sf_latest_2 <- recovered_sf_latest_2 %>%
filter(X != 0) %>%
filter(Y != 0)
recovered_sf_latest_3 <- recovered_sf_latest_2 %>%
select(place, X, Y, value, date) %>%
st_drop_geometry() %>%
spread(date, value)
dates <- colnames(recovered_sf_latest_3)[4:5]
colnames(recovered_sf_latest_3)[4:5] <- c("a", "b")
recovered_sf_latest_3$a[is.na(recovered_sf_latest_3$a)] <- 0
recovered_sf_latest_3$b[is.na(recovered_sf_latest_3$b)] <- 0
recovered_sf_latest_3 <- recovered_sf_latest_3 %>%
mutate(change = b - a)3.4 Plotting
ggplot()+
theme_grey(base_size=8.5)+
theme(axis.title = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size = 13,face="bold"))+
geom_sf(data = countries50, colour = "grey", fill= "white", size=0.2)+
geom_point(data = recovered_sf_latest_3, aes(x = X,
y = Y,
size = change,
colour = change,
alpha= change))+
scale_alpha(range = c(0.1, 1))+
scale_colour_gradientn(colours= c("red1", "red4"))+
scale_size_continuous(range = c(1, 5))+
labs(size= "Recovered cases", alpha = "Recovered cases", colour = "Recovered cases",
title = "Recovered from COVID-19 during the last week")+
guides(alpha = F,
size= F)
recovered_sf_latest_3_change_top10 <- recovered_sf_latest_3 %>%
arrange(-change) %>%
head(n = 10)
ggplot()+
theme_grey(base_size=8.5)+
theme(axis.title = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size = 13,face="bold"),
plot.subtitle = element_text(size=10))+
geom_sf(data = countries50, fill = "white", colour = "grey", size=0.2)+
geom_point(data = recovered_sf_latest_3_change_top10, aes(x = X,
y = Y,
size = change,
alpha = change), colour = "blue")+
scale_size_continuous( range = c(1,6))+
scale_alpha(range = c(0.3, 0.9))+
labs(title="Recovered from COVID-19 during the last week",subtitle = "(Top 10 countries)", size= "Recovered cases", alpha = "Recovered cases", colour = "Recovered cases")
ggplot()+
theme_minimal(base_size=8.5)+
theme(plot.title = element_text(size = 13,face="bold"),
plot.subtitle = element_text(size=10))+
geom_segment(data = recovered_sf_latest_3_change_top10, aes(x= 0,
xend = change,
y = reorder(place, change),
yend = reorder(place, change)),
colour = "black")+
geom_point(data = recovered_sf_latest_3_change_top10, aes(x = change, y = reorder(place, change)),
shape = 21, colour = "white", fill = "firebrick", size= 3)+
labs(title="Recovered from COVID-19 during the last week",subtitle = "(Top 10 countries)",
x = "Recovered cases",
y = "Country")
From graphs we can see that number of recovered cases vary significantly (from 0 to 40000). Europe and North America have most recovered cases, but it would be more meaningful if we examine recovered cases per capita or per cases. The fact that, for example, US and Turkey have the most recovered cases does not means that these countries outperform other countries.
3.5 Population density map for Georgia using classes
ge_adm_0 <- st_read("geo_admbnda_adm0_geostat_20191018.shp")## Reading layer `geo_admbnda_adm0_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm0_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 13 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 40.00834 ymin: 41.05455 xmax: 46.72672 ymax: 43.5865
## geographic CRS: WGS 84ge_adm_1 <- st_read("geo_admbnda_adm1_geostat_20191018.shp")## Reading layer `geo_admbnda_adm1_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm1_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 13 features and 16 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 40.00834 ymin: 41.05455 xmax: 46.72672 ymax: 43.5865
## geographic CRS: WGS 84ge_adm_2 <- st_read("geo_admbnda_adm2_geostat_20191018.shp")## Reading layer `geo_admbnda_adm2_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm2_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 70 features and 19 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 41.54717 ymin: 41.05455 xmax: 46.72672 ymax: 43.25688
## geographic CRS: WGS 84ge_pop<-read_xlsx("geo_admpop.xlsx",sheet=3)## New names:
## * `` -> ...8ge_pop_2 <- ge_pop%>%
select(ADM2_PCODE, POP_Total)
ge_adm_2_ii <- ge_adm_2 %>%
select(ADM2_PCODE, ADM2_EN)
ge_adm_2_pop <- left_join(ge_adm_2_ii, ge_pop_2, by="ADM2_PCODE")
ge_adm_1_TBI <- ge_adm_1 %>%
select(ADM1_PCODE, ADM1_EN) %>%
rename(ADM_PCODE = ADM1_PCODE,
ADM_EN = ADM1_EN) %>%
filter(ADM_EN == "Tbilisi")
ge_adm_1_TBI <- ge_adm_1_TBI %>%
mutate(POP_Total=1171100)
ge_adm_2_pop <- ge_adm_2_pop %>%
rename(ADM_PCODE = ADM2_PCODE,
ADM_EN = ADM2_EN)
ge_adm_2_pop <- rbind(ge_adm_2_pop, ge_adm_1_TBI)ge_adm_2_pop$area <- st_area(ge_adm_2_pop)
ge_adm_2_pop$area <- set_units(ge_adm_2_pop$area, km^2)ge_adm_2_pop <- ge_adm_2_pop %>%
mutate(pop_dens =as.numeric( POP_Total / area))
classes <- classIntervals(ge_adm_2_pop$pop_dens, n = 10, style = "quantile", dig.lab=4)
ge_adm_2_pop <- ge_adm_2_pop %>%
mutate(percent_class = cut(pop_dens, classes$brks, include.lowest = T, dig.lab = 4))ggplot()+
theme_minimal(base_size=8.5)+
theme(axis.title = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size = 13,face="bold"))+
geom_sf(data = ge_adm_0, fill = "grey40", colour = "grey60", size = 0.2)+
geom_sf(data = ge_adm_2_pop, aes(fill = percent_class), colour = "grey60", size= 0.2)+
scale_fill_brewer(type=seq,palette="RdBu",direction=-1) +
labs(title = "Population Density in Georgia",
fill = "Population/Km2",
caption = "data: The Humanitarian Data Exchange \nvisual: D.Kajaia")+
annotation_scale(location = "bl", style = "ticks") +
annotation_north_arrow(location = "tr", width = unit(0.5, "cm"))
I think population density distribution in Georgia depends on some life conditions. For example, we can see that density is the lowest in the northern regions of Georgia where landscape is hilly and density is the highest in big cities where GDP per capita is very high relative to other regions.
4 Part 4
4.1 Libraries
library(rnaturalearth)
library(tidyverse)
library(tmaptools)
library(lubridate)
library(ggspatial)
library(ggthemes)
library(maptools)
library(GISTools)
library(ggplot2)
library(raster)
library(rgeos)
library(tmap)
library(maps)
library(sf)
4.2 Importing data
setwd("C:/Users/User/Documents/OneDrive/Desktop/R/spatial")prec <- raster::getData('worldclim', var ='bio', res = 0.5, lon = 42, lat = 44)
gadm_0 <- raster::getData('GADM', country = 'GEO', level = 0)
gadm_0_sf <- st_as_sf(gadm_0)
prec2 <- mask(prec$bio12_17, gadm_0)
roads <- ne_download(scale = "large", type = "roads", category = "cultural")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\User\AppData\Local\Temp\Rtmp0YCM6x", layer: "ne_10m_roads"
## with 56601 features
## It has 31 fields
## Integer64 fields read as strings: scalerank questionroads_ge <- gIntersection(roads, gadm_0, byid = TRUE, drop_lower_td = TRUE)download.file("http://aasa.ut.ee/tbilisi/osm_georgia_selected.zip", destfile="osm_georgia_selected.zip")
unzip("osm_georgia_selected.zip")
file.remove("osm_georgia_selected.zip")rivers <- st_read("gis_osm_waterways_free_1.shp")## Reading layer `gis_osm_waterways_free_1' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\gis_osm_waterways_free_1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 17123 features and 5 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 39.52382 ymin: 40.93108 xmax: 46.8273 ymax: 45.74497
## geographic CRS: WGS 84rivers_2 <- rivers %>% filter(fclass == "river")
rivers_2 <- st_intersection(rivers_2, gadm_0_sf)## although coordinates are longitude/latitude, st_intersection assumes that they are planar## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriespopulated_places <- ne_download(scale = "large", type = "populated_places", category = "cultural")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\User\AppData\Local\Temp\Rtmp0YCM6x", layer: "ne_10m_populated_places"
## with 7343 features
## It has 119 fields
## Integer64 fields read as strings: wof_id ne_idpopulated_places_sf <- st_as_sf(populated_places)
populated_places_sf <- populated_places_sf %>% filter(ISO_A2 =="GE")4.3 Annual Precipitation in Georgia
plot(gadm_0)
plot(prec2, add= T, col = rainbow(20))
plot(rivers_2, add=T, col="lightblue", lwd = 0.25)
plot(roads_ge, add = T, col = "red", lwd = 0.1)
plot(populated_places_sf, col = "black", add = T)
pointLabel(st_coordinates(populated_places_sf),labels = populated_places_sf$NAME,cex=1)
north.arrow(xb=45.75, yb=43.2, len=0.05, lab="N")
maps::map.scale(x=41.5, y=41.15, relwidth = 0.15,ratio=FALSE, cex=0.9)
title(main ="Annual Precipitation in Georgia", font.main= 1)
4.4 Annual Precipitation in Georgia with ggplot2
re<-as.data.frame(prec2,xy=TRUE,na.rm=TRUE)
str(re)## 'data.frame': 109577 obs. of 3 variables:
## $ x : num 40.2 40.2 40.2 40.3 40.3 ...
## $ y : num 43.6 43.6 43.6 43.6 43.6 ...
## $ bio12_17: int 1265 1223 1200 1215 1226 1223 1279 1319 1259 1285 ...re <- st_as_sf(re, coords = c("x", "y"), crs = 4326)
roads2<-st_as_sf(roads_ge)ggplot()+
theme_grey(base_size=8.5)+
theme(axis.title = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size = 13,face="bold"))+
geom_sf(data = re, aes(color=bio12_17))+
geom_sf(data=roads2,color="red",size=0.1)+
geom_sf(data=rivers_2,color="lightblue")+
geom_sf_label(populated_places_sf,mapping=aes(label=NAME),alpha = 0.6,
size= 3,expand=TRUE,label.size=0.1)+
scale_color_gradientn(colours = rainbow(16))+
labs(color = "mm",
title = "Annual Precipitation in Georgia",
caption = "Author: D. Kajaia")+
annotation_scale(location = "bl", style = "ticks") +
annotation_north_arrow(location = "tr", width = unit(0.5, "cm"))## Warning: Ignoring unknown parameters: expand## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
4.5 Average annual precipitation in Georgian regions
gadm_1_sp <- raster::getData('GADM', country = 'GEO', level = 1)
region_Ave_prec <- raster::extract(prec2, gadm_1_sp, fun = mean, na.rm=TRUE, sp = T)
region_Ave_prec_sf <- st_as_sf(region_Ave_prec)
head(region_Ave_prec_sf)## Simple feature collection with 6 features and 11 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 40.01111 ymin: 41.03851 xmax: 46.72136 ymax: 43.58454
## CRS: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## GID_0 NAME_0 GID_1 NAME_1 VARNAME_1 NL_NAME_1 TYPE_1
## 1 GEO Georgia GEO.1_1 Abkhazia Sokhumi <NA> Avtonomiuri Respublika
## 5 GEO Georgia GEO.2_1 Ajaria Batumi <NA> Avtonomiuri Respublika
## 6 GEO Georgia GEO.3_1 Guria Ozurgeti <NA> Region
## 7 GEO Georgia GEO.4_1 Imereti Kutaisi <NA> Region
## 8 GEO Georgia GEO.5_1 Kakheti Telavi <NA> Region
## 9 GEO Georgia GEO.6_1 Kvemo Kartli Rustavi <NA> Region
## ENGTYPE_1 CC_1 HASC_1 bio12_17 geometry
## 1 Autonomous Republic <NA> GE.AB 1274.5940 MULTIPOLYGON (((41.9095 42....
## 5 Autonomous Republic <NA> GE.AJ 1267.8889 MULTIPOLYGON (((41.93469 41...
## 6 Region <NA> GE.GU 1589.8551 MULTIPOLYGON (((41.72903 42...
## 7 Region <NA> GE.IM 1089.0719 MULTIPOLYGON (((42.67832 41...
## 8 Region <NA> GE.KA 708.0507 MULTIPOLYGON (((45.7786 41....
## 9 Region <NA> GE.KK 623.4208 MULTIPOLYGON (((44.56081 41...ggplot()+
geom_sf(data = region_Ave_prec_sf, aes(fill = bio12_17), col = "grey", size= 0.25)+
geom_sf_label(data = region_Ave_prec_sf, aes(label= round(bio12_17, 0)), alpha = 0.5)+
scale_fill_gradientn(colours = c("#B0B579","#383b23"))+
labs(title = "Average annual precipitation in Georgian regions", fill = "mm")## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
5 Part 5
5.1 Libraries
library(data.table)
library(tidyverse)
library(lubridate)
library(ggspatial)
library(ggthemes)
library(scico)
library(ggmap)
library(knitr)
library(dplyr)
library(sf)
5.2 Importing GPS data
setwd("C:/Users/User/Documents/OneDrive/Desktop/R/spatial")url <- "http://aasa.ut.ee/Rspatial/data/usa_GPS.zip"
download.file(url, "usa_GPS.zip")
unzip("usa_GPS.zip")
file.remove("usa_GPS.zip")gpsData <- read.csv2("gps_us.csv")5.3 Cleaning data
## Rows: 185,808
## Columns: 17
## $ user_id <int> 241, 241, 241, 241, 241, 241, 241, 241, 241, 241, 24...
## $ id <int> 646380266, 646380267, 646380268, 646380269, 64638027...
## $ restart <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2...
## $ counter <dbl> 66718, 66720, 66723, 66730, 66732, 66734, 66736, 667...
## $ time_system <dbl> 1.524345e+12, 1.524345e+12, 1.524345e+12, 1.524345e+...
## $ time_gps <dbl> 1.524345e+12, 1.524345e+12, 1.524345e+12, 1.524345e+...
## $ time_system_ts <chr> "2018-04-22 00:10:39", "2018-04-22 00:10:40", "2018-...
## $ time_gps_ts <chr> "2018-04-22 00:10:36", "2018-04-22 00:10:37", "2018-...
## $ accuracy <dbl> 12, 16, 16, 6, 4, 6, 6, 6, 4, 4, 6, 6, 6, 6, 6, 6, 6...
## $ altitude <dbl> -23.991000, -24.298600, -12.608500, -9.616940, -8.17...
## $ bearing <dbl> 165.957, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ speed <dbl> 0.104187, NA, 0.000000, 0.000000, 0.000000, 0.000000...
## $ header_id <int> 226334354, 226334354, 226334354, 226334360, 22633436...
## $ series_id <int> 1042, 1042, 1042, 1042, 1042, 1042, 1042, 1042, 1042...
## $ year <int> 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018...
## $ X <dbl> -121.3842, -121.3842, -121.3842, -121.3842, -121.384...
## $ Y <dbl> 38.64978, 38.64976, 38.64976, 38.64976, 38.64971, 38...## time_system_ts accuracy altitude bearing speed X Y
## 1 2018-04-22 00:10:39 12 -23.99100 165.957 0.104187 -121.3842 38.64978
## 2 2018-04-22 00:10:40 16 -24.29860 NA NA -121.3842 38.64976
## 3 2018-04-22 00:10:41 16 -12.60850 NA 0.000000 -121.3842 38.64976
## 4 2018-04-22 00:10:42 6 -9.61694 NA 0.000000 -121.3842 38.64976
## 5 2018-04-22 00:10:43 4 -8.17016 NA 0.000000 -121.3843 38.64971
## 6 2018-04-22 00:10:44 6 -5.67983 NA 0.000000 -121.3843 38.649715.4 Importing base map
box <- c(left = min(gpsData2$X)-0.2, bottom = min(gpsData2$Y)-0.2, right = max(gpsData2$X)+0.2, top = max(gpsData2$Y)+0.2)
Ca_map <- get_stamenmap(box, maptype = "terrain", zoom = 9)5.5 Plot of movement path
ggmap(Ca_map)+
theme_classic()+
theme(axis.title = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size = 13,face="bold"))+
geom_point(data = gpsData2, aes(y = Y, x=X), colour ="red", size=0.6, alpha=0.5)+
labs(title ="Movement path on 2018-04-27")+
annotation_north_arrow(location = "tr", width = unit(0.5, "cm"))
6 Part 6
6.1 Libraries
library(readxl)
library(dplyr)
library(knitr)
library(curl)
library(tmap)
library(sf)6.2 Importing data
link<-"http://aasa.ut.ee/tbilisi/data/tripadvisor.xlsx"
curl_download(link,destfile="data.xlsx")data<-read_xlsx("data.xlsx")6.3 Tidying data
glimpse(data)## Rows: 987
## Columns: 4
## $ Title <chr> "Shef's Grandma", "1001 Nights Restaurant", "11 Katkha", "1...
## $ lon <chr> "44.795958900000002", "44.802948100000002", "44.76982009999...
## $ lat <chr> "41.695589300000002", "41.693018000000002", "41.72188049999...
## $ reviews <dbl> 35, 8, 1, 97, 45, 82, 28, 2, 2, 16, 2, 4, 194, 81, 81, 1, 2...data$lon<-as.numeric(data$lon)## Warning: NAs introduced by coerciondata$lat<-as.numeric(data$lat)## Warning: NAs introduced by coerciondata2<-data%>%filter(!is.na(lat))%>%filter(!is.na(lon))
dim(data)==dim(data2)## [1] FALSE TRUE6.4 Converting to sf
data_sf<-st_as_sf(data2, coords = c("lon", "lat"), crs = 4326)
glimpse(data_sf)## Rows: 966
## Columns: 3
## $ Title <chr> "Shef's Grandma", "1001 Nights Restaurant", "11 Katkha", "...
## $ reviews <dbl> 35, 8, 1, 97, 45, 82, 28, 2, 2, 16, 2, 4, 194, 81, 81, 1, ...
## $ geometry <POINT [°]> POINT (44.79596 41.69559), POINT (44.80295 41.69302)...6.5 Map of restaurants in Tbilisi
tmap_mode("view")## tmap mode set to interactive viewingtm_shape(data_sf)+
tm_dots(size = "reviews", col = "red", alpha = 0.5, scale = 1, border.col="black", sizes.legend=c(0.5,2))## Legend for symbol sizes not available in view mode.7 Part 7
7.1 Libraries
library(rnaturalearth)
library(tidyverse)
library(wbstats)
library(dplyr)
library(tmap)
7.2 Choosing indicator
wbstats_indicators <- wb_cachelist$indicators
write.csv(wbstats_indicators, file = "wbstats_indicators.csv")After exporting indicators to CSV files I chose “Renewable electricity (% in total electricity output): Share of electrity generated by renewable power plants in total electricity” with ID - “4.1_SHARE.RE.IN.ELECTRICITY”
data<- wb(indicator ="4.1_SHARE.RE.IN.ELECTRICITY",
mrv = 1, return_wide = TRUE)
glimpse(data)## Rows: 232
## Columns: 5
## $ iso3c <chr> "ABW", "AFG", "AGO", "ALB", "AND", "A...
## $ date <chr> "2015", "2015", "2015", "2015", "2015...
## $ iso2c <chr> "AW", "AF", "AO", "AL", "AD", "AE", "...
## $ country <chr> "Aruba", "Afghanistan", "Angola", "Al...
## $ `4.1_SHARE.RE.IN.ELECTRICITY` <dbl> 14.8561614, 86.0501113, 53.1749283, 1...colnames(data)[5]<-"Value"
glimpse(data)## Rows: 232
## Columns: 5
## $ iso3c <chr> "ABW", "AFG", "AGO", "ALB", "AND", "ARE", "ARG", "ARM", "AS...
## $ date <chr> "2015", "2015", "2015", "2015", "2015", "2015", "2015", "20...
## $ iso2c <chr> "AW", "AF", "AO", "AL", "AD", "AE", "AR", "AM", "AS", "AG",...
## $ country <chr> "Aruba", "Afghanistan", "Angola", "Albania", "Andorra", "Un...
## $ Value <dbl> 14.8561614, 86.0501113, 53.1749283, 100.0000000, 86.1167002...7.3 Downloading map
countries <- ne_download(scale = 50, type = 'countries', category = 'cultural', returnclass = "sf")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\User\AppData\Local\Temp\Rtmp0YCM6x", layer: "ne_50m_admin_0_countries"
## with 241 features
## It has 94 fields
## Integer64 fields read as strings: POP_EST NE_IDcountries2 <- countries %>%
dplyr::select(ADM0_A3, NAME, GDP_MD_EST, POP_EST)7.4 Joining data table with geospatial layer
data_joined <- left_join(countries2, data, by = c("ADM0_A3" = "iso3c"))7.5 Creating plots
tmap_mode("plot")## tmap mode set to plottingtm_shape(data_joined)+
tm_polygons(col = "Value",
style = "kmeans",
n=5,
palette=c('#d7191c','#fdae61','#ffffbf','#a6d96a','#1a9641'),
border.col = "grey50",
lwd = 0.2,
title = "Percentage",
colorNA = "grey80")+
tm_style("natural")+
tm_format("NLD")+
tm_layout(main.title = "The Share of Renewable Electricity",
legend.bg.color = "white",
legend.bg.alpha = 0.6)+
tm_credits("Data: World Bank\nMap: D. Kajaia",
bg.color = "white",
bg.alpha = 0.7)
tmap_mode("view")## tmap mode set to interactive viewingtm_shape(data_joined)+
tm_polygons(col = "Value",
style = "kmeans",
n=5,
palette=c('#d7191c','#fdae61','#ffffbf','#a6d96a','#1a9641'),
border.col = "grey50",
lwd = 0.2,
title = "Percentage",
colorNA = "grey80")+
tm_style("natural")+
tm_format("NLD")+
tm_layout(main.title = "The Share of Renewable Electricity",
legend.bg.color = "white",
legend.bg.alpha = 0.6)+
tm_credits("Data: World Bank\nMap: D. Kajaia",
bg.color = "white",
bg.alpha = 0.7)## Credits not supported in view mode.## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs).
## Need '+proj=longlat +datum=WGS84'